% Figure 4 and 5
axoptions = {'scaled ticks = false',...
             'y tick label style={/pgf/number format/.cd, fixed, fixed zerofill, precision=0, set thousands separator={}}',...
             'x tick label style={/pgf/number format/.cd, precision=1, set thousands separator={}}',...
             'legend style={font=\normalsize}'}; %\scriptsize
% Load data
load ../data/DataforMatlab/IP.csv
load ../data/DataforMatlab/torn.csv
load ../data/DataforMatlab/dec_torn.csv
 
%% Percentile for plotting
run("../subroutines/sub_make_percentile.m")

dec_plot= dec_torn(T,:);
Dec_infN = zeros(N,1);  

% calculations for plotting
    for p = 1:10
        if p == 1
            Dec_infN(1:PC(p,T),1) = dec_plot(p);
        elseif p == 10
            Dec_infN(PC(p-1,T):end,1) = dec_plot(p);
        else        
            Dec_infN(PC(p-1,T):PC(p,T),1) = dec_plot(p);
        end
    end

% To reduce the size of the tikzit figure, select a sample and plot it.
sample_ind = [indmin(T):20:indmax(T),indmax(T)];
OurxaxisUKplot = W(sample_ind,end)*52;
OuryaxisUKplot = log(W(sample_ind,end)./U(sample_ind,end)) - torn(end,1);
OuryaxisUKdes = log(W(sample_ind,end)./U(sample_ind,end)) - Dec_infN(sample_ind);

%% Plotting
figure('name','Figure 5 (a)','NumberTitle','off')
semilogx(OurxaxisUKplot,OuryaxisUKplot*100,'LineWidth', 2);
xlim([min(OurxaxisUKplot),100000])
xticklabels({'10,000','100,000'})
ylabel('Bias (log points)')
xlabel('Income in 2017')
%matlab2tikz('../fig/Bias_MM.tex','extraAxisOptions',axoptions);

% Figure 5 (b)
figure('name','Figure 5 (b)','NumberTitle','off')
semilogx(OurxaxisUKplot, OuryaxisUKdes * 100, 'LineWidth', 2);
xlim([min(OurxaxisUKplot), 100000])
xticklabels({'10,000', '100,000'})
ylabel('Bias (log points)')
xlabel('Income in 2017')
%matlab2tikz('../fig/Bias_Dec.tex', 'extraAxisOptions', axoptions);

% Figure 4 (a)
% To reduce the size of the tikzit figure, select a sample and plot it.
RC = [Wmin;Wq;Wmax] ./ exp(torn') * 52;
OurxaxisUK = [Wmin(:,end);Wq(:,end);Wmax(:,end)] * 52;
OuryaxisUK = [exp(logUmin(:,end)+log(52)); (Uq(:,end) * 52); exp(logUmax(:,end)+log(52))];

figure('name','Figure 4 (a)','NumberTitle','off')
h = loglog(OurxaxisUK, OuryaxisUK, 'LineWidth', 2); hold on
loglog(OurxaxisUK, (RC(:,end)), 'LineWidth', 2);

ax = ancestor(h, 'axes');
ax.XAxis.Exponent = 0;
ax.XAxis.TickLabelFormat = '%.0f';
xlim([min(OurxaxisUK), 100000])
ylim([min(OuryaxisUK), max(OuryaxisUK)])
xticklabels({'10,000', '100,000'})
yticklabels({'1,000', '10,000'})
xlabel('Income in 2017')
ylabel('1974 Income that gives the same utility as in 2017')
legend('Money metric', 'Real consumption', 'Location', 'northwest')
%matlab2tikz('../fig/Money.tex', 'extraAxisOptions', axoptions);

%% save for missing price analysis
save("../data/DataforMatlab/OuryaxisUK_base","OurxaxisUK","OuryaxisUK") % for missing price application
save("../data/DataforMatlab/Uq_base.mat","Uq") % for missing price application

Bx_vec = sum(B_vec([1:9 11:16],:,:),1);
UTMed =quantile(U_vec(:,:,T),0.5);
for t = 1:T
Fbt_U = interpNaN(log(U_vec(:,:,t)),Bx_vec(:,:,t));
bt_U_med(:,t) = Fbt_U(log(UTMed));
end
bt_U_med_base = bt_U_med;
save("../data/DataforMatlab/Compensated_base.mat","bt_U_med_base")

    function F_X = interpNaN(X,Y)
        indexNaN = isnan(X);
        X(indexNaN) = [];% Eliminate NaN
        Y(indexNaN) = [];% Eliminate NaN
        [X, indexUN] = unique(X); 
        F_X = griddedInterpolant(X,Y(indexUN));
    end

